home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_400 / 418_02 / rasmol2 / molecule.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-03-02  |  58.6 KB  |  2,089 lines
  1. /* molecule.c
  2.  * RasMol2 Molecular Graphics
  3.  * Roger Sayle, February 1994
  4.  * Version 2.3
  5.  */
  6. #include "rasmol.h"
  7.  
  8. #ifdef IBMPC
  9. #include <windows.h>
  10. #include <malloc.h>
  11. #endif
  12. #ifndef sun386
  13. #include <stdlib.h>
  14. #endif
  15.  
  16. #include <string.h>
  17. #include <ctype.h>
  18. #include <stdio.h>
  19. #include <math.h>
  20.  
  21. #define MOLECULE
  22. #include "molecule.h"
  23. #include "command.h"
  24. #include "abstree.h"
  25. #include "transfor.h"
  26. #include "render.h"
  27.  
  28.  
  29. #define GroupPool    8
  30. #define HBondPool   32
  31. #define BondPool    32
  32. #define AtomPool    32
  33.  
  34. #define NoLadder     0x00
  35. #define ParaLadder   0x01
  36. #define AntiLadder   0x02
  37.  
  38. #define Cos70Deg     0.34202014332567
  39.  
  40. #define FeatHelix    1
  41. #define FeatSheet    2
  42. #define FeatTurn     3
  43.  
  44. #define SourceNone   0
  45. #define SourcePDB    1
  46. #define SourceCalc   2
  47.  
  48. typedef struct {
  49.         int init, term;
  50.         char chain;
  51.         char type;
  52.         } FeatEntry;
  53.  
  54. #define FeatSize    32
  55. typedef struct _Feature {
  56.     struct _Feature __far *fnext;
  57.         FeatEntry data[FeatSize];
  58.         int count;
  59.         } Feature;
  60.  
  61. typedef struct {
  62.           char src[4];
  63.           char dst[4];
  64.           } ConvTable;
  65.  
  66. #define MAXALCATOM   20
  67. static ConvTable AlcAtomTable[MAXALCATOM] = {
  68.     { { 'C', '3', ' ', ' ' }, { ' ', 'C', '3', ' ' } },  /*  0 */
  69.     { { 'O', '3', ' ', ' ' }, { ' ', 'O', '3', ' ' } },  /*  1 */
  70.     { { 'N', '3', ' ', ' ' }, { ' ', 'N', '3', ' ' } },  /*  2 */
  71.     { { 'S', 'O', '2', ' ' }, { ' ', 'S', '2', ' ' } },  /*  3 */
  72.     { { 'C', 'A', 'R', ' ' }, { ' ', 'C', ' ', ' ' } },  /*  4 */
  73.     { { 'C', '2', ' ', ' ' }, { ' ', 'C', '2', ' ' } },  /*  5 */
  74.     { { 'C', '1', ' ', ' ' }, { ' ', 'C', '1', ' ' } },  /*  6 */
  75.     { { 'O', '2', ' ', ' ' }, { ' ', 'O', '2', ' ' } },  /*  7 */
  76.     { { 'N', '2', ' ', ' ' }, { ' ', 'N', '2', ' ' } },  /*  8 */
  77.     { { 'N', '1', ' ', ' ' }, { ' ', 'N', '1', ' ' } },  /*  9 */
  78.     { { 'N', 'A', 'R', ' ' }, { ' ', 'N', ' ', ' ' } },  /* 10 */
  79.     { { 'N', 'A', 'M', ' ' }, { ' ', 'N', ' ', ' ' } },  /* 11 */
  80.     { { 'N', 'P', 'L', '3' }, { ' ', 'N', '3', ' ' } },  /* 12 */
  81.     { { 'N', '3', '+', ' ' }, { ' ', 'N', '3', ' ' } },  /* 13 */
  82.     { { 'P', '3', ' ', ' ' }, { ' ', 'P', '3', ' ' } },  /* 14 */
  83.     { { 'P', '3', 'D', ' ' }, { ' ', 'P', '3', ' ' } },  /* 15 */
  84.     { { 'P', '4', ' ', ' ' }, { ' ', 'P', '4', ' ' } },  /* 16 */
  85.     { { 'S', '3', ' ', ' ' }, { ' ', 'S', '3', ' ' } },  /* 17 */
  86.     { { 'S', '2', ' ', ' ' }, { ' ', 'S', '2', ' ' } },  /* 18 */
  87.     { { 'S', 'O', ' ', ' ' }, { ' ', 'S', ' ', ' ' } },  /* 19 */
  88.                                  };
  89.  
  90. static Molecule __far *FreeMolecule;
  91. static HBond __far *FreeHBond;
  92. static Chain __far *FreeChain;
  93. static Group __far *FreeGroup;
  94. static Atom __far *FreeAtom;
  95. static Bond __far *FreeBond;
  96.  
  97. static Feature __far *FeatList;
  98. static HBond __far * __far *CurHBond;
  99. static Molecule __far *CurMolecule;
  100. static Chain __far *CurChain;
  101. static Group __far *CurGroup;
  102. static Atom __far *CurAtom;
  103. static Atom __far *ConnectAtom;
  104. static int InfoStrucSource;
  105. static int HMinMaxFlag;
  106. static int MMinMaxFlag;
  107. static int MemSize;
  108.  
  109. static char Record[82];
  110. static FILE *DataFile;
  111.  
  112. /* Macros for commonly used loops */
  113. #define ForEachAtom  for(chain=Database->clist;chain;chain=chain->cnext) \
  114.                      for(group=chain->glist;group;group=group->gnext)    \
  115.                      for(aptr=group->alist;aptr;aptr=aptr->anext)
  116. #define ForEachBond  for(bptr=Database->blist;bptr;bptr=bptr->bnext)
  117.  
  118.  
  119.  
  120. static void FatalDataError(ptr)
  121.     char *ptr;
  122. {
  123.     char buffer[80];
  124.  
  125.     sprintf(buffer,"Database Error: %s!",ptr);
  126.     RasMolFatalExit(buffer);
  127. }
  128.  
  129.  
  130. static void FetchRecord()
  131. {
  132.     register char *ptr;
  133.     register int ch;
  134.  
  135.     ptr = Record + 1;
  136.     if( !feof(DataFile) )
  137.     {   do {
  138.             ch = getc(DataFile);
  139.             if( (ch=='\n') || (ch=='\r') )
  140.             {   if( ptr != Record+1 )
  141.                 {   *ptr = 0;
  142.                     return;
  143.                 }
  144.             } else if( ch==EOF )
  145.             {   *ptr = 0;
  146.                 return;
  147.             } else *ptr++ = ch;
  148.         } while( ptr < Record+80 );
  149.  
  150.         /* skip to the end of the line! */
  151.         do { ch = getc(DataFile);
  152.         } while( (ch!='\n') && (ch!='\r') && (ch!=EOF) );
  153.     }
  154.     *ptr = 0;
  155. }
  156.  
  157.  
  158. static void ExtractString( len, src, dst )
  159.     int len;  char *src, *dst;
  160. {
  161.     register char *ptr;
  162.     register char ch;
  163.     register int i;
  164.  
  165.     ptr = dst;
  166.     for( i=0; i<len; i++ )
  167.     {   if( (ch = *src++) )
  168.         {   if( (*dst++ = ch) != ' ' ) 
  169.                 ptr = dst;
  170.         } else break;
  171.     }
  172.     *ptr = 0;
  173. }
  174.  
  175.  
  176. static Long ReadValue( pos, len )
  177.     int pos, len;
  178. {
  179.     register Long result;
  180.     register char *ptr;
  181.     register char ch;
  182.     register int neg;
  183.  
  184.     result = 0;
  185.     neg = False;
  186.     ptr = Record+pos;
  187.     while( len-- )
  188.     {   ch = *ptr++;
  189.         if( (ch>='0') && (ch<='9') )
  190.         {   result = (10*result)+(ch-'0');
  191.         } else if( ch=='-' )
  192.             neg = True;
  193.     }
  194.     return( neg? -result : result );
  195. }
  196.  
  197.  
  198. void DescribeMolecule()
  199. {
  200.     char buffer[40];
  201.  
  202.     if( CommandActive )
  203.         WriteChar('\n');
  204.     CommandActive=False;
  205.  
  206.     if( *InfoMoleculeName )
  207.     {   WriteString("Molecule name ....... ");
  208.         WriteString(InfoMoleculeName);
  209.         WriteChar('\n');
  210.     }
  211.  
  212.     if( *InfoClassification )
  213.     {   WriteString("Classification ...... ");
  214.         WriteString(InfoClassification);
  215.         WriteChar('\n');
  216.     }
  217.  
  218.     if( Database )
  219.     {   WriteString("Secondary Structure . ");
  220.         if( InfoStrucSource==SourceNone )
  221.         {   WriteString("No Assignment\n");
  222.         } else if( InfoStrucSource==SourcePDB )
  223.         {   WriteString("PDB Data Records\n");
  224.         } else WriteString("Calculated\n");
  225.     }
  226.  
  227.  
  228.     if( *InfoIdentCode )
  229.     {   WriteString("Brookhaven Code ..... ");
  230.         WriteString(InfoIdentCode);
  231.         WriteChar('\n');
  232.     }
  233.  
  234.     if( InfoChainCount>1 )
  235.     {   sprintf(buffer,"Number of Chains .... %d\n",InfoChainCount);
  236.         WriteString(buffer);
  237.     }
  238.  
  239.     sprintf(buffer,"Number of Groups .... %d",MainGroupCount);
  240.     WriteString(buffer);
  241.     if( HetaAtomCount )
  242.     {   sprintf(buffer," (%d)\n",HetaGroupCount);
  243.         WriteString(buffer);
  244.     } else WriteChar('\n');
  245.  
  246.     sprintf(buffer,"Number of Atoms ..... %ld",MainAtomCount);
  247.     WriteString(buffer);
  248.     if( HetaAtomCount )
  249.     {   sprintf(buffer," (%d)\n",HetaAtomCount);
  250.         WriteString(buffer);
  251.     } else WriteChar('\n');
  252.  
  253.     if( InfoBondCount )
  254.     {   sprintf(buffer,"Number of Bonds ..... %ld\n",InfoBondCount);
  255.         WriteString(buffer);
  256.     }
  257.  
  258.     if( InfoSSBondCount != -1 )
  259.     {   WriteString("Number of Bridges ... ");
  260.         sprintf(buffer,"%d\n\n",InfoSSBondCount);
  261.         WriteString(buffer);
  262.     }
  263.  
  264.     if( InfoHBondCount != -1 )
  265.     {   WriteString("Number of H-Bonds ... ");
  266.         sprintf(buffer,"%d\n",InfoHBondCount);
  267.         WriteString(buffer);
  268.     }
  269.  
  270.     if( InfoHelixCount != -1 )
  271.     {   WriteString("Number of Helices ... ");
  272.         sprintf(buffer,"%d\n",InfoHelixCount);
  273.         WriteString(buffer);
  274.  
  275.         WriteString("Number of Ladders ... ");
  276.         sprintf(buffer,"%d\n",InfoLadderCount);
  277.         WriteString(buffer);
  278.  
  279.         WriteString("Number of Turns ..... ");
  280.         sprintf(buffer,"%d\n",InfoTurnCount);
  281.         WriteString(buffer);
  282.     }
  283. }
  284.  
  285.  
  286. static void CreateChain( ident )
  287.     char ident;
  288. {
  289.     register Chain __far *prev;
  290.  
  291.     if( !CurMolecule )
  292.     {   if( !(CurMolecule = FreeMolecule) )
  293.         {   MemSize += sizeof(Molecule);
  294.             CurMolecule = (Molecule __far *)_fmalloc(sizeof(Molecule));
  295.             if( !CurMolecule ) FatalDataError("Memory allocation failed");
  296.         } else FreeMolecule = (void __far*)0;
  297.  
  298.         CurChain = (void __far*)0;
  299.         CurMolecule->slist = (void __far*)0;
  300.         CurMolecule->hlist = (void __far*)0;
  301.         CurMolecule->blist = (void __far*)0;
  302.         Database = CurMolecule;
  303.     }
  304.  
  305.     prev = CurChain;
  306.     if( !(CurChain = FreeChain) )
  307.     {   MemSize += sizeof(Chain);
  308.         CurChain = (Chain __far *)_fmalloc(sizeof(Chain));
  309.         if( !CurChain ) FatalDataError("Memory allocation failed");
  310.     } else FreeChain = FreeChain->cnext;
  311.  
  312.     if( prev )
  313.     {   prev->cnext = CurChain;
  314.     } else CurMolecule->clist = CurChain;
  315.     CurChain->cnext = (void __far*)0;
  316.  
  317.     CurChain->ident = ident;
  318.     CurChain->glist = (void __far*)0;
  319.     CurChain->blist = (void __far*)0;
  320.     ConnectAtom = (void __far*)0;
  321.     CurGroup = (void __far*)0;
  322.     InfoChainCount++;
  323. }
  324.  
  325.  
  326. static Atom __far *CreateAtom()
  327. {
  328.     register Atom __far *ptr;
  329.     register int i;
  330.  
  331.     if( !(ptr = FreeAtom) )
  332.     {   MemSize += AtomPool*sizeof(Atom);
  333.         ptr = (Atom __far *)_fmalloc( AtomPool*sizeof(Atom) );
  334.         if( !ptr ) FatalDataError("Memory allocation failed");
  335.         for( i=1; i<AtomPool; i++ )
  336.         {   ptr->anext = FreeAtom;
  337.             FreeAtom = ptr++;
  338.         } 
  339.     } else FreeAtom = ptr->anext;
  340.  
  341.     if( CurAtom )
  342.     {   CurAtom->anext = ptr;
  343.     } else CurGroup->alist = ptr;
  344.     ptr->anext = (void __far*)0;
  345.     CurAtom = ptr;
  346.  
  347.     SelectCount++;
  348.     ptr->flag = SelectFlag;
  349.     ptr->radius = 375;
  350.     ptr->altl = ' ';
  351.     ptr->mbox = 0;
  352.     ptr->col = 0;
  353.  
  354.     return( ptr );
  355. }
  356.  
  357.  
  358. static void ProcessAtom( ptr )
  359.     Atom __far *ptr;
  360. {
  361.     if( GetElemIdent(ptr) == HandH )
  362.         ptr->flag |= HydrogenFlag;
  363.     if( !(ptr->flag&(HydrogenFlag|HeteroFlag)) )
  364.         ptr->flag |= NormAtomFlag;
  365.  
  366. #ifdef INVERT
  367.     ptr->yorg = -ptr->yorg;
  368. #endif
  369.  
  370.     if( HMinMaxFlag || MMinMaxFlag )
  371.     {   if( ptr->xorg < MinX ) 
  372.         {   MinX = ptr->xorg;
  373.         } else if( ptr->xorg > MaxX ) 
  374.             MaxX = ptr->xorg;
  375.  
  376.         if( ptr->yorg < MinY ) 
  377.         {   MinY = ptr->yorg;
  378.         } else if( ptr->yorg > MaxY ) 
  379.             MaxY = ptr->yorg;
  380.  
  381.         if( ptr->zorg < MinZ ) 
  382.         {   MinZ = ptr->zorg;
  383.         } else if( ptr->zorg > MaxZ ) 
  384.             MaxZ = ptr->zorg;
  385.     } else 
  386.     {   MinX = MaxX = ptr->xorg;
  387.         MinY = MaxY = ptr->yorg;
  388.         MinZ = MaxZ = ptr->zorg;
  389.     }
  390.             
  391.     if( ptr->flag & HeteroFlag )
  392.     {   if( HMinMaxFlag )
  393.         {   if( ptr->temp < MinHetaTemp ) 
  394.             {   MinHetaTemp = ptr->temp;
  395.             } else if( ptr->temp > MaxHetaTemp ) 
  396.                 MaxHetaTemp = ptr->temp;
  397.         } else MinHetaTemp = MaxHetaTemp = ptr->temp;
  398.         HMinMaxFlag = True;
  399.         HetaAtomCount++;
  400.     } else
  401.     {   if( MMinMaxFlag )
  402.         {   if( ptr->temp < MinMainTemp ) 
  403.             {   MinMainTemp = ptr->temp;
  404.             } else if( ptr->temp > MaxMainTemp ) 
  405.                 MaxMainTemp = ptr->temp;
  406.         } else MinMainTemp = MaxMainTemp = ptr->temp;
  407.         MMinMaxFlag = True;
  408.         MainAtomCount++;
  409.     }
  410. }
  411.  
  412.  
  413. static int FindResNo( ptr )
  414.     char *ptr;
  415. {
  416.     register int refno;
  417.  
  418.     for( refno=0; refno<ResNo; refno++ )
  419.         if( !strncmp(Residue[refno],ptr,3) )
  420.             return( refno );
  421.  
  422.     if( !strncmp(ptr,"CSH",3) ||
  423.         !strncmp(ptr,"CYH",3) ||
  424.         !strncmp(ptr,"CSM",3) )
  425.     {   return(17);  /* cystine */
  426.     } else if( !strncmp(ptr,"WAT",3) ||
  427.                !strncmp(ptr,"H20",3) )
  428.     /* check HHO, OHH, 0H2 and SOL? */
  429.     {   return(31);
  430.     } else if( !strncmp(ptr,"D20",3) )
  431.     {   return(32);  /* Heavy water (DOD) */
  432.     } else if( !strncmp(ptr,"SUL",3) )
  433.     {   return(33);  /* Sulphate (SO4) */
  434.     } else if( !strncmp(ptr,"CPR",3) )
  435.     {   return(11);  /* cis-proline */
  436.     } else if( !strncmp(ptr,"TRY",3) )
  437.         return(15);  /* tryptophan */
  438.     
  439.  
  440.     if( ResNo++ == MAXRES )
  441.         FatalDataError("Too many new residues");
  442.     Residue[refno][0] = *ptr++;
  443.     Residue[refno][1] = *ptr++;
  444.     Residue[refno][2] = *ptr;
  445.     return( refno );
  446. }
  447.  
  448.  
  449. static void ProcessPDBColourMask()
  450. {
  451.     register MaskDesc *ptr;
  452.     register char *mask;
  453.     register int i;
  454.  
  455.     if( MaskCount==MAXMASK )
  456.         FatalDataError("Too many COLOR records in file");
  457.     ptr = &UserMask[MaskCount];
  458.     mask = ptr->mask;
  459.  
  460.  
  461.     ptr->flags = 0;
  462.     for( i=7; i<12; i++ )
  463.         if( (*mask++ = Record[i]) != '#' )
  464.             ptr->flags |= SerNoFlag;
  465.  
  466.     for( i=13; i<21; i++ )
  467.         *mask++ = Record[i];
  468.     *mask++ = Record[22];
  469.  
  470.     for( i=23; i<27; i++ )
  471.         if( (*mask++ = Record[i]) != '#' )
  472.             ptr->flags |= ResNoFlag;
  473.     *mask++ = Record[27];
  474.  
  475.     ptr->r = (int)(ReadValue(31,8)>>2) + 5;
  476.     ptr->g = (int)(ReadValue(39,8)>>2) + 5;
  477.     ptr->b = (int)(ReadValue(47,8)>>2) + 5;
  478.     ptr->radius = (short)(5*ReadValue(55,6))>>1;
  479.     MaskCount++;
  480. }
  481.  
  482.  
  483. static Bond __far *ProcessBond( src, dst, flag )
  484.     Atom __far *src, __far *dst;
  485.     int flag;
  486. {
  487.     register Bond __far *ptr;
  488.     register int i;
  489.  
  490.     if( !(ptr = FreeBond) )
  491.     {   MemSize += BondPool*sizeof(Bond);
  492.         ptr = (Bond __far *)_fmalloc( BondPool*sizeof(Bond) );
  493.         if( !ptr ) FatalDataError("Memory allocation failed");
  494.         for( i=1; i<BondPool; i++ )
  495.         {   ptr->bnext = FreeBond;
  496.             FreeBond = ptr++;
  497.         } 
  498.     } else FreeBond = ptr->bnext;
  499.     
  500.     ptr->flag = flag | SelectFlag;
  501.     ptr->srcatom = src;
  502.     ptr->dstatom = dst;
  503.     ptr->radius = 0;
  504.     ptr->col = 0;
  505.  
  506.     return( ptr );
  507. }
  508.  
  509.  
  510.  
  511. static void ProcessPDBGroup( heta, serno )
  512.     int heta, serno;
  513. {
  514.     register Group __far *ptr;
  515.     register int i;
  516.  
  517.     if( !CurChain || (CurChain->ident!=Record[22]) )
  518.         CreateChain( Record[22] );
  519.  
  520.     if( !(ptr = FreeGroup) )
  521.     {   MemSize += GroupPool*sizeof(Group);
  522.         ptr = (Group __far *)_fmalloc( GroupPool*sizeof(Group) );
  523.         if( !ptr ) FatalDataError("Memory allocation failed");
  524.         for( i=1; i<GroupPool; i++ )
  525.         {   ptr->gnext = FreeGroup;
  526.             FreeGroup = ptr++;
  527.         } 
  528.     } else FreeGroup = ptr->gnext;
  529.     
  530.     if( CurGroup )
  531.     {   CurGroup->gnext = ptr;
  532.     } else CurChain->glist = ptr;
  533.     ptr->gnext = (void __far*)0;
  534.     CurGroup = ptr;
  535.  
  536.     CurAtom = (void __far*)0;
  537.     ptr->refno = FindResNo( &Record[18] );
  538.     ptr->alist = (void __far*)0;
  539.     ptr->serno = serno;
  540.     ptr->flag = 0;
  541.     ptr->col = 0;
  542.  
  543.     /* Solvents should be hetero! */
  544.     if( IsSolvent(ptr->refno) )
  545.         heta = True;
  546.  
  547.     if( heta )
  548.     {   HetaGroupCount++;
  549.         if( HMinMaxFlag )
  550.         {   if( serno > MaxHetaRes ) 
  551.             {   MaxHetaRes = serno;
  552.             } else if( serno < MinHetaRes )
  553.                 MinHetaRes = serno;
  554.         } else MinHetaRes = MaxHetaRes = serno;
  555.     } else 
  556.     {   MainGroupCount++;
  557.         if( MMinMaxFlag )
  558.         {   if( serno > MaxMainRes )
  559.             {   MaxMainRes = serno;
  560.             } else if( serno < MinHetaRes )
  561.                 MinHetaRes = serno;
  562.         } else MinMainRes = MaxMainRes = serno; 
  563.     }
  564. }
  565.  
  566.  
  567. static void ProcessPDBAtom( heta )
  568.     int heta;
  569. {
  570.     register Bond __far *bptr;
  571.     register Atom __far *ptr;
  572.     register Long dx,dy,dz;
  573.     register int serno;
  574.     register int refno;
  575.     register int i;
  576.  
  577.     /* Ignore Pseudo Atoms!! */
  578.     if( (Record[13]==' ') && (Record[14]=='Q') )
  579.         return; 
  580.  
  581.     serno = (int)ReadValue(23,4);
  582.     if( !CurGroup || (CurGroup->serno!=serno) 
  583.         || (CurChain->ident!=Record[22]) )
  584.         ProcessPDBGroup( heta, serno );
  585.  
  586.     /* Solvents should be hetero! */
  587.     if( IsSolvent(CurGroup->refno) )
  588.         heta = True;
  589.  
  590.  
  591.     ptr = CreateAtom();
  592.     for( refno=0; refno<ElemNo; refno++ )
  593.         if( !strncmp(ElemDesc[refno],Record+13,4) )
  594.             break;
  595.  
  596.     if( refno==ElemNo )
  597.     {   if( ElemNo++ == MAXELEM )
  598.             FatalDataError("Too many new atom types");
  599.  
  600.         for( i=0; i<4; i++ )
  601.             ElemDesc[refno][i] = Record[i+13];
  602.     }
  603.  
  604.     ptr->refno = refno;
  605.     ptr->altl = Record[17];
  606.     ptr->serno = (int)ReadValue(7,5);
  607.     ptr->temp = (int)ReadValue(61,6);
  608.  
  609.     ptr->xorg =  ReadValue(31,8)/4;
  610.     ptr->yorg =  ReadValue(39,8)/4;
  611.     ptr->zorg = -ReadValue(47,8)/4;
  612.  
  613.     if( heta ) 
  614.         ptr->flag |= HeteroFlag;
  615.  
  616.     ProcessAtom( ptr );
  617.     if( IsAlphaCarbon(refno) && IsAmino(CurGroup->refno) )
  618.     {   if( ConnectAtom )
  619.         {   dx = ConnectAtom->xorg - ptr->xorg;
  620.             dy = ConnectAtom->yorg - ptr->yorg;
  621.             dz = ConnectAtom->zorg - ptr->zorg;
  622.  
  623.             /* Break backbone if CA-CA > 7.00A */
  624.             if( dx*dx+dy*dy+dz*dz < (Long)1750*1750 )
  625.             {   bptr = ProcessBond(ConnectAtom,ptr,NormBondFlag);
  626.                 bptr->bnext = CurChain->blist;
  627.                 CurChain->blist = bptr;
  628.             } else ptr->flag |= BreakFlag;
  629.         }
  630.         ConnectAtom = ptr;
  631.     } else if( IsSugarPhosphate(refno) && IsNucleo(CurGroup->refno) )
  632.     {   if( ConnectAtom )
  633.         {   bptr = ProcessBond(ConnectAtom,ptr,NormBondFlag);
  634.             bptr->bnext = CurChain->blist;
  635.             CurChain->blist = bptr;
  636.         }
  637.         ConnectAtom = ptr;
  638.     }
  639. }
  640.  
  641.  
  642. static FeatEntry __far *AllocFeature()
  643. {
  644.     register Feature __far *ptr;
  645.  
  646.     if( !FeatList || (FeatList->count==FeatSize) )
  647.     {   ptr = (Feature __far*)_fmalloc(sizeof(Feature));
  648.         if( !ptr ) FatalDataError("Memory allocation failed");
  649.         ptr->fnext = FeatList;
  650.         ptr->count = 0;
  651.         FeatList = ptr;
  652.     } else ptr = FeatList;
  653.  
  654.     return( &(ptr->data[ptr->count++]) );
  655. }
  656.  
  657.  
  658. #if defined(__STDC__) || defined(IBMPC)
  659. static void UpdateFeature( FeatEntry __far*, int );
  660. #endif
  661.  
  662. static void UpdateFeature( ptr, mask )
  663.     FeatEntry __far *ptr;  int mask;
  664. {
  665.     register Chain __far *chain;
  666.     register Group __far *group;
  667.  
  668.     for( chain=Database->clist; chain; chain=chain->cnext )
  669.         if( chain->ident == ptr->chain )
  670.         {   group=chain->glist;
  671.             while( group && (group->serno<ptr->init) )
  672.                 group = group->gnext;
  673.  
  674.             while( group && (group->serno<=ptr->term) )
  675.             {   group->flag |= mask;
  676.                 group = group->gnext;
  677.             }
  678.             return;
  679.         }
  680. }
  681.  
  682.  
  683. static void ProcessPDBFeatures()
  684. {
  685.     register Feature __far *next;
  686.     register Feature __far *ptr;
  687.     register int i;
  688.  
  689.     InfoTurnCount = 0;
  690.     InfoHelixCount = 0;
  691.     InfoLadderCount = 0;
  692.     InfoStrucSource = SourcePDB;
  693.  
  694.     for( ptr=FeatList; ptr; ptr=next )
  695.     {    if( Database )
  696.              for( i=0; i<ptr->count; i++ )
  697.                  if( ptr->data[i].type==FeatHelix )
  698.                  {   UpdateFeature( &ptr->data[i], HelixFlag );
  699.                      InfoHelixCount++;
  700.                  } else if( ptr->data[i].type==FeatSheet )
  701.                  {   UpdateFeature( &ptr->data[i], SheetFlag );
  702.                      InfoLadderCount++;
  703.                  } else /* FeatTurn */
  704.                  {   UpdateFeature( &ptr->data[i], TurnFlag );
  705.                      InfoTurnCount++;
  706.                  }
  707.  
  708.          /* Deallocate Memory */
  709.          next = ptr->fnext;
  710.          _ffree( ptr );
  711.     }
  712. }
  713.  
  714.  
  715. int LoadPDBMolecule( fp )
  716.     FILE *fp;
  717. {
  718.     register FeatEntry __far *ptr;
  719.     register char *src, *dst;
  720.     register int model;
  721.  
  722.     model = True;
  723.     FeatList = (void __far*)0;
  724.     DataFile = fp;
  725.  
  726.  
  727.     do {   
  728.         FetchRecord();
  729.  
  730.         if( !strncmp("ATOM",Record+1,4) )
  731.         {   if( model ) ProcessPDBAtom(False);
  732.         } else if( !strncmp("HETA",Record+1,4) )
  733.         {   if( model ) ProcessPDBAtom(True);
  734.         } else if( !strncmp("SHEE",Record+1,4) )
  735.         {   if( !model ) continue;
  736.         
  737.             /* Remaining SHEET record fields   */
  738.             /* 39-40 .... Strand Parallelism   */
  739.             /* 33 ....... Same Chain as 22?    */
  740.             ptr = AllocFeature();
  741.             ptr->type = FeatSheet;
  742.             ptr->chain = Record[22];
  743.             ptr->init = (int)ReadValue(23,4);
  744.             ptr->term = (int)ReadValue(34,4);
  745.        
  746.         } else if( !strncmp("HELI",Record+1,4) )
  747.         {   if( !model ) continue;
  748.         
  749.             /* Remaining HELIX record fields   */
  750.             /* 39-40 .... Helix Classification */
  751.             /* 32 ....... Same Chain as 20?    */
  752.             ptr = AllocFeature();
  753.             ptr->type = FeatHelix;
  754.             ptr->chain = Record[20];
  755.             ptr->init = (int)ReadValue(22,4);
  756.             ptr->term = (int)ReadValue(34,4);
  757.  
  758.         } else if( !strncmp("TURN",Record+1,4) )
  759.         {   if( !model ) continue;
  760.         
  761.             ptr = AllocFeature();
  762.             ptr->type = FeatTurn;
  763.             ptr->chain = Record[20];
  764.             ptr->init = (int)ReadValue(21,4);
  765.             ptr->term = (int)ReadValue(32,4);
  766.             
  767.         } else if( !strncmp("COLO",Record+1,4) )
  768.         {   ProcessPDBColourMask();
  769.         } else if( !strncmp("TER ",Record+1,4) )
  770.         {   ConnectAtom = (void __far*)0;
  771.         
  772.         } else if( !strncmp("HEAD",Record+1,4) )
  773.         {   ExtractString(40,Record+11,InfoClassification);
  774.             ExtractString( 4,Record+63,InfoIdentCode);
  775.             AdviseUpdate(AdvClass);
  776.             AdviseUpdate(AdvIdent);
  777.             
  778.         } else if( !strncmp("COMP",Record+1,4) )
  779.         {   if( Record[10]==' ' )  /* First COMPND record */
  780.             {   ExtractString(60,Record+11,InfoMoleculeName);
  781.                 AdviseUpdate(AdvName);
  782.             }
  783.         } else if( !strncmp("CRYS",Record+1,4) )
  784.         {   dst = InfoSpaceGroup;
  785.             for( src=Record+56; *src && src<Record+67; src++ )
  786.                 if( *src!=' ' ) *dst++ = *src;
  787.             *dst = 0;
  788.         } else if( !strncmp("ENDM",Record+1,4) )
  789.         {   model = False;
  790.         } else if( !strncmp("END ",Record+1,4) ) 
  791.             break;
  792.     } while( !feof(DataFile) );
  793.  
  794.     if( Database )
  795.         strcpy(InfoFileName,DataFileName);
  796.     if( FeatList ) ProcessPDBFeatures();
  797.     return( True );
  798. }
  799.  
  800.  
  801. int LoadXYZMolecule( fp )
  802.     FILE *fp;
  803. {
  804.     return( False );
  805. }
  806.  
  807.  
  808. static int FindAlchemyResNo()
  809. {
  810.     register char *ptr;
  811.     register int i, j;
  812.     char name[4];
  813.  
  814.     if( Record[8]==' ' )
  815.     {   name[0] = name[2] = name[3] = ' ';
  816.         if( islower(Record[7]) )
  817.         {   name[1] = toupper(Record[7]); 
  818.         } else name[1] = Record[7];
  819.         ptr = name;
  820.     } else
  821.     {   ptr = Record+7;
  822.         for( i=0; i<4; i++ )
  823.         if( islower(ptr[i]) )
  824.                 ptr[i] = toupper(ptr[i]);
  825.  
  826.         for( i=0; i<MAXALCATOM; i++ )
  827.             if( !strncmp(AlcAtomTable[i].src,ptr,4) )
  828.             {   ptr = AlcAtomTable[i].dst;
  829.                 break;
  830.             }
  831.     }
  832.  
  833.     for( j=0; j<ElemNo; j++ )
  834.         if( !strncmp(ElemDesc[j],ptr,4) )
  835.             return( j );
  836.  
  837.     if( ElemNo++ == MAXELEM )
  838.         FatalDataError("Too many new atom types");
  839.     for( i=0; i<4; i++ ) ElemDesc[j][i] = ptr[i];
  840.     return( j );
  841. }
  842.  
  843.  
  844. int LoadAlchemyMolecule( fp )
  845.     FILE *fp;
  846. {
  847.     register Atom __far *ptr;
  848.     register int atoms, bonds;
  849.     register int i;
  850.  
  851.     DataFile = fp;
  852.     FetchRecord();
  853.     atoms = (int)ReadValue(1,5);
  854.     bonds = (int)ReadValue(14,5);
  855.     ExtractString(38,Record+42,InfoMoleculeName);
  856.  
  857.     if( atoms )
  858.     {   CreateChain( ' ' );
  859.         if( !(CurGroup = FreeGroup) )
  860.         {   MemSize += sizeof(Group);
  861.             CurGroup = (Group __far *)_fmalloc(sizeof(Group));
  862.             if( !CurGroup ) FatalDataError("Memory allocation failed");
  863.         } else FreeGroup = CurGroup->gnext;
  864.         strcpy(InfoFileName,DataFileName);
  865.  
  866.         CurAtom = (void __far*)0;
  867.         CurChain->glist = CurGroup;
  868.         CurGroup->refno = FindResNo( "MOL" );
  869.         CurGroup->gnext = (void __far*)0;
  870.         CurGroup->alist = (void __far*)0;
  871.         CurGroup->serno = 1;
  872.         CurGroup->flag = 0;
  873.         CurGroup->col = 0;
  874.         
  875.         MinMainRes = MaxMainRes = 1;
  876.         MinHetaRes = MaxHetaRes = 0;
  877.         MainGroupCount = 1;
  878.  
  879.         for( i=0; i<atoms; i++ )
  880.         {   FetchRecord();
  881.             ptr = CreateAtom();
  882.  
  883.             ptr->refno = FindAlchemyResNo();
  884.             ptr->temp = (int)ReadValue(41,8);
  885.             /* ptr->serno = (int)ReadValue(1,5); */
  886.             ptr->serno = i;
  887.  
  888.             ptr->xorg = ReadValue(13,7)/4;
  889.             ptr->yorg = ReadValue(22,7)/4;
  890.             ptr->zorg = ReadValue(31,7)/4;
  891.             ProcessAtom( ptr );
  892.         }
  893.  
  894.         for( i=0; i<bonds; i++ )
  895.         {   /* InfoBondCount++; */
  896.             FetchRecord();
  897.         }
  898.     }
  899.     return( True );
  900. }
  901.  
  902.  
  903.  
  904. int SavePDBMolecule( filename )
  905.     char *filename;
  906. {
  907.     register Real x, y, z;
  908.     register float xpos, ypos, zpos;
  909.     register Group __far *prev;
  910.     register Chain __far *chain;
  911.     register Group __far *group;
  912.     register Atom __far *aptr;
  913.     register char *ptr;
  914.     register int count;
  915.     register char ch;
  916.     register int i;
  917.  
  918.     if( !Database )
  919.         return( False );
  920.  
  921.     DataFile = fopen( filename, "w" );
  922.     if( !DataFile )
  923.     {   if( CommandActive )
  924.             WriteChar('\n');
  925.         WriteString("Error: Unable to create file!\n\n");
  926.         CommandActive=False;
  927.         return( False );
  928.     }
  929.  
  930.     if( *InfoClassification || *InfoIdentCode )
  931.     {   fputs("HEADER    ",DataFile);
  932.  
  933.         ptr = InfoClassification;
  934.         for( i=11; i<=50; i++ )
  935.             putc( (*ptr ? *ptr++ : ' '), DataFile );
  936.         fprintf(DataFile,"13-JUL-93   %.4s\n",InfoIdentCode);
  937.     }
  938.  
  939.     if( *InfoMoleculeName )
  940.         fprintf(DataFile,"COMPND    %.60s\n",InfoMoleculeName);
  941.  
  942.     prev = (void __far*)0;
  943.  
  944.     count = 1;
  945.     ForEachAtom
  946.         if( aptr->flag&SelectFlag )
  947.         {   if( prev && (chain->ident!=ch) )
  948.                 fprintf( DataFile, "TER   %5d      %.3s %c%4d \n", 
  949.                          count++, Residue[prev->refno], ch, prev->serno);
  950.  
  951.             if( aptr->flag&HeteroFlag )
  952.             {      fputs("HETATM",DataFile);
  953.             } else fputs("ATOM  ",DataFile);
  954.             fprintf( DataFile, "%5d %.4s %.3s %c%4d    ",
  955.                      count++, ElemDesc[aptr->refno], Residue[group->refno],
  956.                      chain->ident, group->serno );
  957.  
  958.             x = aptr->xorg/250.0;
  959.             y = aptr->yorg/250.0;
  960.             z = aptr->zorg/250.0;
  961.  
  962.             /* Apply Current Viewpoint Rotation Matrix */
  963.             xpos = (float)(x*RotX[0] + y*RotX[1] + z*RotX[2]);
  964.             ypos = (float)(x*RotY[0] + y*RotY[1] + z*RotY[2]);
  965.             zpos = (float)(x*RotZ[0] + y*RotZ[1] + z*RotZ[2]);
  966.  
  967. #ifdef INVERT
  968.             fprintf(DataFile,"%8.3f%8.3f%8.3f",xpos,-ypos,-zpos);
  969. #else
  970.             fprintf(DataFile,"%8.3f%8.3f%8.3f",xpos,ypos,-zpos);
  971. #endif
  972.             fprintf(DataFile,"  1.00%6.2f\n",aptr->temp/100.0);
  973.             
  974.             ch = chain->ident;
  975.             prev = group;
  976.         }
  977.  
  978.     if( prev )
  979.         fprintf( DataFile, "TER   %5d      %.3s %c%4d \n", 
  980.                  count, Residue[prev->refno], ch, prev->serno);
  981.  
  982.     fputs("END   \n",DataFile);
  983.     fclose( DataFile );
  984.     return( True );
  985. }
  986.  
  987.  
  988. int SaveXYZMolecule( filename )
  989.     char *filename;
  990. {
  991.     return( True );
  992. }
  993.  
  994.  
  995. int SaveAlchemyMolecule( filename )
  996.     char *filename;
  997. {
  998.     register Real x, y, z;
  999.     register float xpos, ypos, zpos;
  1000.     register Chain __far *chain;
  1001.     register Group __far *group;
  1002.     register Atom __far *aptr;
  1003.     register Bond __far *bptr;
  1004.     register char *ptr;
  1005.     register int atomno;
  1006.     register int bondno;
  1007.     register int num;
  1008.  
  1009.     if( !Database )
  1010.         return( False );
  1011.  
  1012.     DataFile = fopen( filename, "w" );
  1013.     if( !DataFile )
  1014.     {   if( CommandActive )
  1015.             WriteChar('\n');
  1016.         WriteString("Error: Unable to create file!\n\n");
  1017.         CommandActive=False;
  1018.         return( False );
  1019.     }
  1020.  
  1021.     atomno = 0;
  1022.     ForEachAtom
  1023.         if( aptr->flag & SelectFlag )
  1024.         {   aptr->mbox = 0;
  1025.             atomno++;
  1026.         }
  1027.  
  1028.     bondno = 0;
  1029.     ForEachBond
  1030.         if( (bptr->srcatom->flag&bptr->dstatom->flag) & SelectFlag )
  1031.         {   if( bptr->flag&AromBondFlag )
  1032.             {   bptr->srcatom->mbox = -1;
  1033.                 bptr->dstatom->mbox = -1;
  1034.             } else if( !(bptr->flag&HydrBondFlag) )
  1035.             {   num = (bptr->flag&DoubBondFlag)? 2 : 1;
  1036.                 if( bptr->srcatom->mbox>0 ) 
  1037.                     bptr->srcatom->mbox += num;
  1038.                 if( bptr->dstatom->mbox>0 ) 
  1039.                     bptr->dstatom->mbox += num;
  1040.             }
  1041.             bondno++;
  1042.         }
  1043.  
  1044.     fprintf(DataFile,"%5d ATOMS, %5d BONDS, ",atomno,bondno);
  1045.     fprintf(DataFile,"    0 CHARGES, %s\n", InfoMoleculeName );
  1046.  
  1047.     atomno = 1;
  1048.     ForEachAtom
  1049.         if( aptr->flag & SelectFlag )
  1050.         {   aptr->mbox = atomno;
  1051.             fprintf(DataFile,"%5d ",atomno++);
  1052.  
  1053.             switch( GetElemIdent(aptr) )
  1054.             {   case( HandC ):   if( aptr->mbox == -1 )
  1055.                                  {   ptr = "CAR ";
  1056.                                  } else if( aptr->mbox == 1 )
  1057.                                  {   ptr = "C3  ";
  1058.                                  } else ptr = "C2  ";
  1059.                                  fputs( ptr, DataFile );
  1060.                                  break;
  1061.  
  1062.                 case( HandN ):   if( aptr->mbox == -1 )
  1063.                                  {   ptr = "NAR ";
  1064.                                  } else ptr = "N2  ";
  1065.                                  fputs( ptr, DataFile );
  1066.                                  break;
  1067.  
  1068.                 case( HandO ):   if( aptr->mbox == 2 )
  1069.                                  {   ptr = "O2  ";
  1070.                                  } else ptr = "O3  ";
  1071.                                  fputs( ptr, DataFile );
  1072.                                  break;
  1073.  
  1074.                 case( HandH ):   fputs( "H   ", DataFile );  break;
  1075.  
  1076.                 default:  ptr = ElemDesc[aptr->refno];
  1077.               if( *ptr==' ' )
  1078.                           {   fprintf(DataFile,"%.3s ",ptr+1);
  1079.                           } else fprintf(DataFile,"%.4s",ptr);
  1080.             }
  1081.  
  1082.             x = aptr->xorg/250.0;
  1083.             y = aptr->yorg/250.0;
  1084.             z = aptr->zorg/250.0;
  1085.  
  1086.             /* Apply Current Viewpoint Rotation Matrix */
  1087.             xpos = (float)(x*RotX[0] + y*RotX[1] + z*RotX[2]);
  1088.             ypos = (float)(x*RotY[0] + y*RotY[1] + z*RotY[2]);
  1089.             zpos = (float)(x*RotZ[0] + y*RotZ[1] + z*RotZ[2]);
  1090.  
  1091. #ifdef INVERT
  1092.             fprintf(DataFile,"  %8.4f %8.4f %8.4f",xpos,-ypos,-zpos);
  1093. #else
  1094.             fprintf(DataFile,"  %8.4f %8.4f %8.4f",xpos,ypos,-zpos);
  1095. #endif
  1096.             fprintf(DataFile,"    %7.4f\n",aptr->temp/1000.0);
  1097.         }
  1098.  
  1099.     bondno = 1;
  1100.     ForEachBond
  1101.         if( (bptr->srcatom->flag&bptr->dstatom->flag) & SelectFlag )
  1102.         {   fprintf(DataFile,"%5d %5d %5d  ", bondno++,
  1103.             bptr->srcatom->mbox, bptr->dstatom->mbox );
  1104.             if( bptr->flag & AromBondFlag )
  1105.             {   ptr = "AROMATIC\n";
  1106.             } else if( bptr->flag & DoubBondFlag )
  1107.             {   ptr = "DOUBLE\n";
  1108.             } else ptr = "SINGLE\n";
  1109.             fputs( ptr, DataFile );
  1110.         }
  1111.  
  1112.     ForEachAtom
  1113.         if( aptr->flag & SelectFlag )
  1114.             aptr->mbox = 0;
  1115.     fclose( DataFile );
  1116.     return( True );
  1117. }
  1118.  
  1119.  
  1120. static void TestBonded( sptr, dptr )
  1121.     Atom __far *sptr, __far *dptr;
  1122. {
  1123.     register Bond __far *bptr;
  1124.     register Long dx, dy, dz;
  1125.     register Long max, min;
  1126.     register Long dist;
  1127.     register int flag;
  1128.  
  1129.  
  1130.     if( ((sptr->flag)&HydrogenFlag) || ((dptr->flag)&HydrogenFlag) )
  1131.     {    flag = HydrBondFlag;
  1132.          max = MaxHBondDist;
  1133.          min = MinHBondDist;
  1134.     } else
  1135.     {    flag = NormBondFlag;
  1136.          max = MaxBondDist;
  1137.          min = MinBondDist;
  1138.     }
  1139.  
  1140.     dx = sptr->xorg-dptr->xorg;   if( (dist=dx*dx)>max ) return;
  1141.     dy = sptr->yorg-dptr->yorg;   if( (dist+=dy*dy)>max ) return;
  1142.     dz = sptr->zorg-dptr->zorg;   if( (dist+=dz*dz)>max ) return;
  1143.     if( dist>=min )  
  1144.     {   /* Reset Non-bonded flags! */
  1145.         sptr->flag &= ~NonBondFlag;
  1146.         dptr->flag &= ~NonBondFlag;
  1147.  
  1148.         bptr = ProcessBond(sptr,dptr,flag);
  1149.         bptr->bnext = CurMolecule->blist;
  1150.         CurMolecule->blist = bptr;
  1151.         InfoBondCount++;
  1152.     }
  1153. }
  1154.  
  1155.  
  1156. static void ReclaimHBonds( ptr )
  1157.     HBond __far *ptr;
  1158. {
  1159.     register HBond __far *temp;
  1160.  
  1161.     if( (temp = ptr) )
  1162.     {   while( temp->hnext )
  1163.             temp=temp->hnext;
  1164.         temp->hnext = FreeHBond;
  1165.         FreeHBond = ptr;
  1166.     }
  1167. }
  1168.  
  1169.  
  1170. static void ReclaimBonds( ptr )
  1171.     Bond __far *ptr;
  1172. {
  1173.     register Bond __far *temp;
  1174.  
  1175.     if( (temp = ptr) )
  1176.     {   while( temp->bnext )
  1177.             temp=temp->bnext;
  1178.         temp->bnext = FreeBond;
  1179.         FreeBond = ptr;
  1180.     }
  1181. }
  1182.  
  1183.  
  1184. void CreateMoleculeBonds( info )
  1185.     int info;
  1186. {
  1187.     register int i, x, y, z;
  1188.     register Long tx, ty, tz;
  1189.     register Long mx, my, mz; 
  1190.     register Long dx, dy, dz;
  1191.     register int lx, ly, lz, hx, hy, hz;
  1192.     register Atom __far *aptr, __far *dptr;
  1193.     register Chain __far *chain;
  1194.     register Group __far *group;
  1195.     char buffer[40];
  1196.  
  1197.  
  1198.     if( !Database ) 
  1199.         return;
  1200.  
  1201.     dx = (MaxX-MinX)+1;
  1202.     dy = (MaxY-MinY)+1;
  1203.     dz = (MaxZ-MinZ)+1;
  1204.  
  1205.     ReclaimBonds( CurMolecule->blist );
  1206.     ResetVoxelData();
  1207.  
  1208.     ForEachAtom
  1209.     {   /* Initially non-bonded! */
  1210.         aptr->flag |= NonBondFlag;
  1211.  
  1212.         mx = aptr->xorg-MinX;
  1213.         my = aptr->yorg-MinY;
  1214.         mz = aptr->zorg-MinZ;
  1215.  
  1216.         tx = mx-AbsMaxBondDist;  
  1217.         ty = my-AbsMaxBondDist;  
  1218.         tz = mz-AbsMaxBondDist;  
  1219.  
  1220.         lx = (tx>0)? (int)((VOXORDER*tx)/dx) : 0;
  1221.         ly = (ty>0)? (int)((VOXORDER*ty)/dy) : 0;
  1222.         lz = (tz>0)? (int)((VOXORDER*tz)/dz) : 0;
  1223.  
  1224.         tx = mx+AbsMaxBondDist;  
  1225.         ty = my+AbsMaxBondDist;  
  1226.         tz = mz+AbsMaxBondDist;
  1227.  
  1228.         hx = (tx<dx)? (int)((VOXORDER*tx)/dx) : VOXORDER-1;
  1229.         hy = (ty<dy)? (int)((VOXORDER*ty)/dy) : VOXORDER-1;
  1230.         hz = (tz<dz)? (int)((VOXORDER*tz)/dz) : VOXORDER-1;
  1231.         
  1232.         for( x=lx; x<=hx; x++ )
  1233.         {   i = VOXORDER2*x + VOXORDER*ly;
  1234.             for( y=ly; y<=hy; y++ )
  1235.             {   for( z=lz; z<=hz; z++ )
  1236.                     if( (dptr = (Atom __far*)HashTable[i+z]) )
  1237.                         do { TestBonded(aptr,dptr);
  1238.                         } while( (dptr = dptr->next) );
  1239.                 i += VOXORDER;
  1240.             }
  1241.         }
  1242.                 
  1243.         x = (int)((VOXORDER*mx)/dx);
  1244.         y = (int)((VOXORDER*my)/dy);
  1245.         z = (int)((VOXORDER*mz)/dz);
  1246.  
  1247.         i = VOXORDER2*x + VOXORDER*y + z;
  1248.         aptr->next = (Atom __far*)HashTable[i];
  1249.         HashTable[i] = (void __far*)aptr;
  1250.     }
  1251.     VoxelsClean = False;
  1252.  
  1253.     if( info )
  1254.     {   if( CommandActive )
  1255.             WriteChar('\n');
  1256.         CommandActive=False;
  1257.         sprintf(buffer,"Number of Bonds ..... %ld\n\n",InfoBondCount);
  1258.         WriteString(buffer);
  1259.     }
  1260. }
  1261.  
  1262.  
  1263. Atom __far *FindGroupAtom( group, n )
  1264.     Group __far *group;  Byte n;
  1265. {
  1266.     register Atom __far *ptr;
  1267.  
  1268.     ptr = group->alist;
  1269.     while( ptr )
  1270.     {   if( ptr->refno == n )
  1271.             return( ptr );
  1272.         ptr = ptr->anext;
  1273.     }
  1274.     return( (Atom __far*)0 );
  1275. }
  1276.  
  1277.  
  1278. void TestDisulphideBridge( grp1, grp2, cys1 )
  1279.     Group __far *grp1, __far *grp2;  
  1280.     Atom __far *cys1;
  1281. {
  1282.     register HBond __far *ptr;
  1283.     register Atom __far *cys2;
  1284.     register int dx, dy, dz;
  1285.     register Long max,dist;
  1286.  
  1287.     if( !(cys2=FindGroupAtom(grp2,20)) )
  1288.         return;
  1289.  
  1290.     max = (Long)750*750;
  1291.     dx = (int)(cys1->xorg-cys2->xorg);   if( (dist=(Long)dx*dx)>max ) return;
  1292.     dy = (int)(cys1->yorg-cys2->yorg);   if( (dist+=(Long)dy*dy)>max ) return;
  1293.     dz = (int)(cys1->zorg-cys2->zorg);   if( (dist+=(Long)dz*dz)>max ) return;
  1294.  
  1295.     if( !(ptr = FreeHBond) )
  1296.     {   MemSize += sizeof(HBond);
  1297.         ptr = (HBond __far*)_fmalloc(sizeof(HBond));
  1298.         if( !ptr ) FatalDataError("Memory allocation failed");
  1299.     } else FreeHBond = ptr->hnext;
  1300.  
  1301.     ptr->dst = cys1;
  1302.     if( !(ptr->dstCA=FindGroupAtom(grp1,1)) )
  1303.         ptr->dstCA = cys1;
  1304.  
  1305.     ptr->src = cys2;
  1306.     if( !(ptr->srcCA=FindGroupAtom(grp2,1)) )
  1307.         ptr->srcCA = cys2;
  1308.  
  1309.     ptr->offset = 0;
  1310.     ptr->energy = 0;
  1311.     ptr->flag = 0;
  1312.     ptr->col = 0;
  1313.  
  1314.     ptr->hnext = CurMolecule->slist;
  1315.     CurMolecule->slist = ptr;
  1316.  
  1317.     grp1->flag |= CystineFlag;
  1318.     grp2->flag |= CystineFlag;
  1319.     InfoSSBondCount++;
  1320. }
  1321.  
  1322.  
  1323. void FindDisulphideBridges()
  1324. {
  1325.     register Chain __far *chn1;
  1326.     register Chain __far *chn2;
  1327.     register Group __far *grp1;
  1328.     register Group __far *grp2;
  1329.     register Atom __far *cys;
  1330.     char buffer[40];
  1331.  
  1332.     if( !Database ) return;
  1333.     ReclaimHBonds( CurMolecule->slist );
  1334.     InfoSSBondCount = 0;
  1335.  
  1336.     for(chn1=Database->clist;chn1;chn1=chn1->cnext)
  1337.         for(grp1=chn1->glist;grp1;grp1=grp1->gnext)
  1338.             if( IsCysteine(grp1->refno) && (cys=FindGroupAtom(grp1,20)) )
  1339.             {   for(grp2=grp1->gnext;grp2;grp2=grp2->gnext)
  1340.                     if( IsCysteine(grp2->refno) )
  1341.                         TestDisulphideBridge(grp1,grp2,cys);
  1342.  
  1343.                 for(chn2=chn1->cnext;chn2;chn2=chn2->cnext)
  1344.                     for(grp2=chn2->glist;grp2;grp2=grp2->gnext)
  1345.                         if( IsCysteine(grp2->refno) )
  1346.                             TestDisulphideBridge(grp1,grp2,cys);
  1347.             }
  1348.  
  1349.     if( CommandActive )
  1350.         WriteChar('\n');
  1351.     CommandActive=False;
  1352.     
  1353.     sprintf(buffer,"Number of Bridges ... %d\n\n",InfoSSBondCount);
  1354.     WriteString(buffer);
  1355. }
  1356.  
  1357. #if defined(__STDC__) || defined(IBMPC)
  1358. static void CreateHydrogenBond( Atom __far*, Atom __far*,
  1359.                                 Atom __far*, Atom __far*, int, int );
  1360. #endif
  1361.  
  1362. static void CreateHydrogenBond( srcCA, dstCA, src, dst, energy, offset )
  1363.     Atom __far *srcCA, __far *dstCA;
  1364.     Atom __far *src, __far *dst;
  1365.     int energy, offset;
  1366. {
  1367.     register HBond __far *ptr;
  1368.     register int i;
  1369.  
  1370.     if( !(ptr = FreeHBond) )
  1371.     {   MemSize += HBondPool*sizeof(HBond);
  1372.         ptr = (HBond __far *)_fmalloc( HBondPool*sizeof(HBond) );
  1373.         if( !ptr ) FatalDataError("Memory allocation failed");
  1374.         for( i=1; i<HBondPool; i++ )
  1375.         {   ptr->hnext = FreeHBond;
  1376.             FreeHBond = ptr++;
  1377.         } 
  1378.     } else FreeHBond = ptr->hnext;
  1379.  
  1380.     if( (offset>=-128) && (offset<127) )
  1381.     {   ptr->offset = (Char)offset;
  1382.     } else ptr->offset = 0;
  1383.  
  1384.     ptr->src = src;
  1385.     ptr->dst = dst;
  1386.     ptr->srcCA = srcCA;
  1387.     ptr->dstCA = dstCA;
  1388.     ptr->energy = energy;
  1389.     ptr->flag = 0;
  1390.     ptr->col = 0;
  1391.  
  1392.     *CurHBond = ptr;
  1393.     ptr->hnext = (void __far*)0;
  1394.     CurHBond = &ptr->hnext;
  1395.     InfoHBondCount++;
  1396. }
  1397.  
  1398. /* Coupling constant for Electrostatic Energy   */
  1399. /* QConst = -332 * 0.42 * 0.2 * 1000.0 [*250.0] */
  1400. #define QConst (-6972000.0)
  1401. #define MaxHDist ((Long)2250*2250)
  1402. #define MinHDist ((Long)125*125)
  1403.  
  1404.  
  1405. /* Protein Donor Atom Coordinates */
  1406. static int hxorg,hyorg,hzorg;
  1407. static int nxorg,nyorg,nzorg;
  1408. static Atom __far *best1CA;
  1409. static Atom __far *best2CA;
  1410. static Atom __far *best1;
  1411. static Atom __far *best2;
  1412. static Atom __far *optr;
  1413. static int res1,res2;
  1414. static int off1,off2;
  1415.  
  1416.  
  1417. static int CalculateBondEnergy( grp )
  1418.     Group __far *grp;
  1419. {
  1420.     register double dho,dhc;
  1421.     register double dnc,dno;
  1422.  
  1423.     register Atom __far *cptr;
  1424.     register Long dx,dy,dz;
  1425.     register Long dist;
  1426.     register int result;
  1427.  
  1428.     if( !(cptr=FindGroupAtom(grp,2)) )  return(0);
  1429.     if( !(optr=FindGroupAtom(grp,3)) )  return(0);
  1430.  
  1431.     dx = hxorg-optr->xorg;  
  1432.     dy = hyorg-optr->yorg;  
  1433.     dz = hzorg-optr->zorg;
  1434.     dist = dx*dx+dy*dy+dz*dz;
  1435.     if( dist < MinHDist ) 
  1436.         return( -9900 );
  1437.     dho = sqrt((double)dist);
  1438.  
  1439.     dx = hxorg-cptr->xorg;  
  1440.     dy = hyorg-cptr->yorg;  
  1441.     dz = hzorg-cptr->zorg;
  1442.     dist = dx*dx+dy*dy+dz*dz;
  1443.     if( dist < MinHDist ) 
  1444.         return( -9900 );
  1445.     dhc = sqrt((double)dist);
  1446.  
  1447.     dx = nxorg-cptr->xorg;  
  1448.     dy = nyorg-cptr->yorg;  
  1449.     dz = nzorg-cptr->zorg;
  1450.     dist = dx*dx+dy*dy+dz*dz;
  1451.     if( dist < MinHDist ) 
  1452.         return( -9900 );
  1453.     dnc = sqrt((double)dist);
  1454.  
  1455.     dx = nxorg-optr->xorg;  
  1456.     dy = nyorg-optr->yorg;  
  1457.     dz = nzorg-optr->zorg;
  1458.     dist = dx*dx+dy*dy+dz*dz;
  1459.     if( dist < MinHDist ) 
  1460.         return( -9900 );
  1461.     dno = sqrt((double)dist);
  1462.  
  1463.     result = (int)(QConst/dho - QConst/dhc + QConst/dnc - QConst/dno);
  1464.  
  1465.     if( result<-9900 ) 
  1466.     {   return( -9900 );
  1467.     } else if( result>-500 ) 
  1468.         return( 0 );
  1469.  
  1470.     return( result );
  1471.  
  1472. }
  1473.  
  1474.  
  1475. void CalcHydrogenBonds()
  1476. {
  1477.     register int energy, offset, refno;
  1478.     register Chain __far *chn1, __far *chn2;
  1479.     register Group __far *grp1, __far *grp2;
  1480.     register Group __far *best;
  1481.     register Atom __far *ca1;
  1482.     register Atom __far *ca2;
  1483.     register Atom __far *pc1;
  1484.     register Atom __far *po1;
  1485.     register Atom __far *n1;
  1486.     register Long max,dist;
  1487.     register int pos1,pos2;
  1488.     register int dx,dy,dz;
  1489.     register double dco;
  1490.     char buffer[40];
  1491.  
  1492.     if( !Database ) return;
  1493.     ReclaimHBonds( CurMolecule->hlist );
  1494.     CurMolecule->hlist = (void __far*)0;
  1495.     CurHBond = &CurMolecule->hlist;
  1496.     InfoHBondCount = 0;
  1497.  
  1498.     if( MainAtomCount > 10000 )
  1499.     {   if( CommandActive )
  1500.             WriteChar('\n');
  1501.         WriteString("Please wait... ");
  1502.         CommandActive=True;
  1503.     }
  1504.  
  1505.     for(chn1=Database->clist;chn1;chn1=chn1->cnext)
  1506.     {   if( !chn1->glist ) continue;
  1507.  
  1508.         if( IsProtein(chn1->glist->refno) )
  1509.         {   pc1 = po1 = (void __far*)0;
  1510.             pos1 = 0;
  1511.             for(grp1=chn1->glist;grp1;grp1=grp1->gnext)
  1512.             {   pos1++;
  1513.                 if( pc1 && po1 )
  1514.                 {   dx = (int)(pc1->xorg - po1->xorg);
  1515.                     dy = (int)(pc1->yorg - po1->yorg);
  1516.                     dz = (int)(pc1->zorg - po1->zorg);
  1517.                 } else
  1518.                 {   pc1 = FindGroupAtom(grp1,2);
  1519.                     po1 = FindGroupAtom(grp1,3);
  1520.                     continue;
  1521.                 }
  1522.  
  1523.                 pc1 = FindGroupAtom(grp1,2);
  1524.                 po1 = FindGroupAtom(grp1,3);
  1525.  
  1526.                 if( !IsAmino(grp1->refno) || IsProline(grp1->refno) )
  1527.                     continue;
  1528.  
  1529.                 if( !(ca1=FindGroupAtom(grp1,1)) ) continue;
  1530.                 if( !(n1=FindGroupAtom(grp1,0)) )  continue;
  1531.  
  1532.                 dist = (Long)dx*dx + (Long)dy*dy + (Long)dz*dz;
  1533.                 dco = sqrt( (double)dist )/250.0;
  1534.  
  1535.                 nxorg = (int)n1->xorg;   hxorg = nxorg + (int)(dx/dco);
  1536.                 nyorg = (int)n1->yorg;   hyorg = nyorg + (int)(dy/dco);
  1537.                 nzorg = (int)n1->zorg;   hzorg = nzorg + (int)(dz/dco);
  1538.                 res1 = res2 = 0;
  1539.  
  1540.                 for(chn2=Database->clist;chn2;chn2=chn2->cnext)
  1541.                 {   /* Only consider non-empty peptide chains! */
  1542.                     if( !chn2->glist || !IsProtein(chn2->glist->refno) )
  1543.                         continue;
  1544.  
  1545.                     pos2 = 0;
  1546.                     for(grp2=chn2->glist;grp2;grp2=grp2->gnext)
  1547.                     {   pos2++;
  1548.                         if( (grp2==grp1) || (grp2->gnext==grp1) )
  1549.                             continue;
  1550.  
  1551.                         if( !IsAmino(grp2->refno) ) 
  1552.                             continue;
  1553.                         if( !(ca2=FindGroupAtom(grp2,1)) ) 
  1554.                             continue;
  1555.  
  1556.                         dx = (int)(ca1->xorg-ca2->xorg);
  1557.                         if( (dist=(Long)dx*dx) > MaxHDist )
  1558.                             continue;
  1559.  
  1560.                         dy = (int)(ca1->yorg-ca2->yorg);
  1561.                         if( (dist+=(Long)dy*dy) > MaxHDist )
  1562.                             continue;
  1563.  
  1564.                         dz = (int)(ca1->zorg-ca2->zorg);
  1565.                         if( (dist+=(Long)dz*dz) > MaxHDist )
  1566.                             continue;
  1567.  
  1568.                         if( (energy = CalculateBondEnergy(grp2)) )
  1569.                         {   if( chn1 == chn2 )
  1570.                             {   offset = pos1 - pos2;
  1571.                             } else offset = 0;
  1572.  
  1573.                             if( energy<res1 )
  1574.                             {   best2CA = best1CA;  best1CA = ca2;
  1575.                                 best2 = best1;      best1 = optr;
  1576.                                 res2 = res1;        res1 = energy;
  1577.                                 off2 = off1;        off1 = offset;
  1578.                             } else if( energy<res2 )
  1579.                             {   best2CA = ca2;
  1580.                                 best2 = optr;
  1581.                                 res2 = energy;
  1582.                                 off2 = offset;
  1583.                             }
  1584.                         }
  1585.                     }  /* grp2 */
  1586.                 }      /* chn2 */
  1587.  
  1588.                 if( res1 ) 
  1589.                 {   if( res2 ) 
  1590.                         CreateHydrogenBond(ca1,best2CA,n1,best2,res2,off2);
  1591.                     CreateHydrogenBond(ca1,best1CA,n1,best1,res1,off1);
  1592.                 }
  1593.             }
  1594.         } else if( IsNucleo(chn1->glist->refno) )
  1595.             for(grp1=chn1->glist;grp1;grp1=grp1->gnext)
  1596.             {   if( !IsPurine(grp1->refno) ) continue;
  1597.                 /* Find N1 of Purine Group */
  1598.                 if( !(n1=FindGroupAtom(grp1,21)) )
  1599.                     continue;
  1600.  
  1601.                 /* Maximum N1-N3 distance 5A */
  1602.                 refno = NucleicCompl(grp1->refno);
  1603.                 max = (Long)1250*1250;
  1604.                 best = (void __far*)0;
  1605.  
  1606.                 for(chn2=Database->clist;chn2;chn2=chn2->cnext)
  1607.                 {   /* Only consider non-empty peptide chains! */
  1608.                     if( (chn1==chn2) || !chn2->glist || 
  1609.                         !IsNucleo(chn2->glist->refno) )
  1610.                         continue;
  1611.  
  1612.                     for(grp2=chn2->glist;grp2;grp2=grp2->gnext)
  1613.                         if( grp2->refno == (Byte)refno )
  1614.                         {   /* Find N3 of Pyramidine Group */
  1615.                             if( !(ca1=FindGroupAtom(grp2,23)) )
  1616.                                 continue;
  1617.  
  1618.                             dx = (int)(ca1->xorg - n1->xorg);
  1619.                             if( (dist=(Long)dx*dx) >= max ) 
  1620.                                 continue;
  1621.  
  1622.                             dy = (int)(ca1->yorg - n1->yorg);
  1623.                             if( (dist+=(Long)dy*dy) >= max ) 
  1624.                                 continue;
  1625.  
  1626.                             dz = (int)(ca1->zorg - n1->zorg);
  1627.                             if( (dist+=(Long)dz*dz) >= max )
  1628.                                 continue;
  1629.  
  1630.                             best1 = ca1;
  1631.                             best = grp2;
  1632.                             max = dist;
  1633.                         }
  1634.                 }
  1635.  
  1636.                 if( best )
  1637.                 {   /* Find the sugar phosphorous atoms */
  1638.                     ca1 = FindGroupAtom( grp1, 7 );
  1639.                     ca2 = FindGroupAtom( best, 7 );
  1640.  
  1641.                     CreateHydrogenBond( ca1, ca2, n1, best1, 0, 0 );
  1642.                     if( IsGuanine(grp1->refno) )
  1643.                     {   /* Guanine-Cytosine */
  1644.                         if( (ca1=FindGroupAtom(grp1,22)) &&    /* G.N2 */
  1645.                             (ca2=FindGroupAtom(best,26)) )     /* C.O2 */
  1646.                             CreateHydrogenBond( (void __far*)0, (void __far*)0,
  1647.                                                 ca1, ca2, 0, 0 );
  1648.  
  1649.                         if( (ca1=FindGroupAtom(grp1,28)) &&    /* G.O6 */
  1650.                             (ca2=FindGroupAtom(best,24)) )     /* C.N4 */
  1651.                             CreateHydrogenBond( (void __far*)0, (void __far*)0,
  1652.                                                 ca1, ca2, 0, 0 );
  1653.  
  1654.                     } else /* Adenine-Thymine */
  1655.                         if( (ca1=FindGroupAtom(grp1,25)) &&    /* A.N6 */
  1656.                             (ca2=FindGroupAtom(best,27)) )     /* T.O4 */
  1657.                             CreateHydrogenBond( (void __far*)0, (void __far*)0,
  1658.                                                 ca1, ca2, 0, 0 );
  1659.                 }
  1660.             }
  1661.     }
  1662.  
  1663.     if( CommandActive )
  1664.         WriteChar('\n');
  1665.     CommandActive=False;
  1666.     
  1667.     sprintf(buffer,"Number of H-Bonds ... %d\n",InfoHBondCount);
  1668.     WriteString(buffer);
  1669. }
  1670.  
  1671. #if defined(__STDC__) || defined(IBMPC)
  1672. static int IsHBonded( Atom __far*, Atom __far*, HBond __far* );
  1673. static void TestLadder( Chain __far* );
  1674. #endif
  1675.  
  1676. static int IsHBonded( src, dst, ptr )
  1677.     Atom __far *src, __far *dst;
  1678.     HBond __far *ptr;
  1679. {
  1680.     while( ptr && (ptr->srcCA==src) )
  1681.         if( ptr->dstCA == dst )
  1682.         {   return( True );
  1683.         } else ptr=ptr->hnext;
  1684.     return( False );
  1685. }
  1686.  
  1687.  
  1688. static void FindAlphaHelix( pitch, flag )
  1689.     int pitch, flag;
  1690. {
  1691.     register HBond __far *hbond;
  1692.     register Chain __far *chain;
  1693.     register Group __far *group;
  1694.     register Group __far *first;
  1695.     register Group __far *ptr;
  1696.     register Atom __far *srcCA;
  1697.     register Atom __far *dstCA;
  1698.     register int res,dist,prev;
  1699.  
  1700.     /* Protein chains only! */
  1701.     hbond = Database->hlist;
  1702.     for( chain=Database->clist; chain; chain=chain->cnext )
  1703.     if( (first=chain->glist) && IsProtein(first->refno) )
  1704.     {   prev = False; dist = 0;
  1705.         for( group=chain->glist; hbond && group; group=group->gnext )
  1706.         {   if( IsAmino(group->refno) && (srcCA=FindGroupAtom(group,1)) )
  1707.             {   if( dist==pitch )
  1708.                 {   res = False;
  1709.                     dstCA=FindGroupAtom(first,1);
  1710.  
  1711.                     while( hbond && hbond->srcCA == srcCA )
  1712.                     {   if( hbond->dstCA==dstCA ) res=True;
  1713.                         hbond = hbond->hnext;
  1714.                     }
  1715.  
  1716.                     if( res )
  1717.                     {   if( prev )
  1718.                         {   if( !(first->flag&HelixFlag) ) 
  1719.                                 InfoHelixCount++;
  1720.  
  1721.                             ptr = first;
  1722.                             do {
  1723.                                 ptr->flag |= flag;
  1724.                                 ptr = ptr->gnext;
  1725.                             } while( ptr != group );
  1726.                         } else prev = True;
  1727.                     } else prev = False;
  1728.                 } else while( hbond && hbond->srcCA==srcCA )
  1729.                     hbond = hbond->hnext;
  1730.             } else prev = False;
  1731.  
  1732.             if( group->flag&HelixFlag )
  1733.             {   first = group; prev = False; dist = 1;
  1734.             } else if( dist==pitch )
  1735.             {   first = first->gnext;
  1736.             } else dist++;
  1737.         }
  1738.     } else if( first && IsNucleo(first->refno) )
  1739.         while( hbond && !IsAminoBackbone(hbond->src->refno) )
  1740.             hbond = hbond->hnext;
  1741. }
  1742.  
  1743.  
  1744. static Atom __far *cprevi, __far *ccurri, __far *cnexti;
  1745. static HBond __far *hcurri, __far *hnexti;
  1746. static Group __far *curri, __far *nexti;
  1747.  
  1748.  
  1749.  
  1750. static void TestLadder( chain )
  1751.     Chain __far *chain;
  1752. {
  1753.     register Atom __far *cprevj, __far *ccurrj, __far *cnextj;
  1754.     register HBond __far *hcurrj, __far *hnextj;
  1755.     register Group __far *currj, __far *nextj;
  1756.     register int count, result, found;
  1757.  
  1758.     /* Already part of atleast one ladder */
  1759.     found = curri->flag & SheetFlag;
  1760.     nextj = nexti->gnext;
  1761.  
  1762.     hnextj = hnexti;
  1763.     while( hnextj && hnextj->srcCA==cnexti )
  1764.         hnextj = hnextj->hnext;
  1765.  
  1766.     while( True )
  1767.     {   if( nextj )
  1768.             if( IsProtein(chain->glist->refno) )
  1769.             {   count = 1;
  1770.                 do {
  1771.                     cnextj = FindGroupAtom(nextj,1);
  1772.                     if( count == 3 )
  1773.                     {   if( IsHBonded(cnexti,ccurrj,hnexti) &&
  1774.                             IsHBonded(ccurrj,cprevi,hcurrj) )
  1775.                         {   result = ParaLadder;
  1776.                         } else if( IsHBonded(cnextj,ccurri,hnextj) &&
  1777.                                    IsHBonded(ccurri,cprevj,hcurri) )
  1778.                         {   result = ParaLadder;
  1779.                         } else if( IsHBonded(cnexti,cprevj,hnexti) &&
  1780.                                    IsHBonded(cnextj,cprevi,hnextj) )
  1781.                         {   result = AntiLadder;
  1782.                         } else if( IsHBonded(ccurrj,ccurri,hcurrj) &&
  1783.                                    IsHBonded(ccurri,ccurrj,hcurri) )
  1784.                         {   result = AntiLadder;
  1785.                         } else result = NoLadder;
  1786.  
  1787.                         if( result )
  1788.                         {   curri->flag |= SheetFlag;
  1789.                             currj->flag |= SheetFlag;
  1790.                             if( found ) return;
  1791.                             found = True;
  1792.                         }
  1793.                     } else count++;
  1794.  
  1795.                     cprevj = ccurrj; ccurrj = cnextj; 
  1796.                     currj = nextj;   hcurrj = hnextj;
  1797.  
  1798.                     while( hnextj && hnextj->srcCA==cnextj )
  1799.                         hnextj = hnextj->hnext;
  1800.                 } while( (nextj = nextj->gnext) );
  1801.  
  1802.             } else if( IsNucleo(chain->glist->refno) )
  1803.                 while( hnextj && !IsAminoBackbone(hnextj->src->refno) )
  1804.                     hnextj = hnextj->hnext;
  1805.  
  1806.         if( (chain = chain->cnext) ) 
  1807.         {   nextj = chain->glist;
  1808.         } else return;
  1809.     }
  1810. }
  1811.  
  1812.  
  1813. static void FindBetaSheets()
  1814. {
  1815.     register Chain __far *chain;
  1816.     register int ladder;
  1817.     register int count;
  1818.  
  1819.     hnexti = Database->hlist;
  1820.     for( chain=Database->clist; chain; chain=chain->cnext )
  1821.         if( (nexti = chain->glist) )
  1822.             if( IsProtein(nexti->refno) )
  1823.             {   count = 1;
  1824.                 ladder = False;
  1825.                 do {
  1826.                     cnexti = FindGroupAtom(nexti,1);
  1827.  
  1828.                     if( count == 3 )
  1829.                     {   TestLadder( chain );
  1830.                         if( curri->flag & SheetFlag )
  1831.                         {   if( ladder )
  1832.                             {   InfoLadderCount++;
  1833.                             } else ladder = True;
  1834.                         } else ladder = False;
  1835.                     } else count++;
  1836.  
  1837.                     cprevi = ccurri; ccurri = cnexti; 
  1838.                     curri = nexti;   hcurri = hnexti;
  1839.                     while( hnexti && hnexti->srcCA==cnexti )
  1840.                         hnexti = hnexti->hnext;
  1841.                 } while( (nexti = nexti->gnext) );
  1842.  
  1843.             } else if( IsNucleo(nexti->refno) )
  1844.                 while( hnexti && !IsAminoBackbone(hnexti->src->refno) )
  1845.                     hnexti = hnexti->hnext;
  1846. }
  1847.  
  1848.  
  1849. static void FindTurnStructure()
  1850. {
  1851.     static Atom __far *aptr[5];
  1852.     register Chain __far *chain;
  1853.     register Group __far *group;
  1854.     register Group __far *prev;
  1855.     register Long ux,uy,uz,mu;
  1856.     register Long vx,vy,vz,mv;
  1857.     register int i,found,len;
  1858.     register Real CosKappa;
  1859.  
  1860.     for( chain=Database->clist; chain; chain=chain->cnext )
  1861.         if( chain->glist && IsProtein(chain->glist->refno) )
  1862.         {   len = 0;  found = False;
  1863.             for( group=chain->glist; group; group=group->gnext )
  1864.             {    if( group->flag & BreakFlag )
  1865.                  {   found = False;
  1866.                      len = 0;
  1867.                  } else if( len==5 )
  1868.                  {   for( i=0; i<4; i++ )
  1869.                          aptr[i] = aptr[i+1];
  1870.                      len = 4;
  1871.                  } else if( len==2 )
  1872.                      prev = group;
  1873.  
  1874.                  aptr[len++] = FindGroupAtom(group,1);
  1875.                  if( len==5 ) 
  1876.                  {   if( !(prev->flag&(HelixFlag|SheetFlag)) &&
  1877.                          aptr[0] && aptr[2] && aptr[4] )
  1878.                      {   ux = aptr[2]->xorg - aptr[0]->xorg;
  1879.                          uy = aptr[2]->yorg - aptr[0]->yorg;
  1880.                          uz = aptr[2]->zorg - aptr[0]->zorg;
  1881.  
  1882.                          vx = aptr[4]->xorg - aptr[2]->xorg;
  1883.                          vy = aptr[4]->yorg - aptr[2]->yorg;
  1884.                          vz = aptr[4]->zorg - aptr[2]->zorg;
  1885.  
  1886.                          mu = ux*ux + uy*uy + uz*uz;
  1887.                          mv = vx*vx + vz*vz + vy*vy;
  1888.                          if( mu && mv )
  1889.                          {   CosKappa = (Real)(ux*vx + uy*vy + uz*vz);
  1890.                              CosKappa /= sqrt( (Real)mu*mv );
  1891.                              if( CosKappa<Cos70Deg )
  1892.                              {   if( !found )
  1893.                                      InfoTurnCount++;
  1894.                                  prev->flag |= TurnFlag;
  1895.                              }
  1896.                          }
  1897.                      }
  1898.                      found = prev->flag&TurnFlag;
  1899.                      prev = prev->gnext;
  1900.                  } /* len==5 */
  1901.             }
  1902.         }
  1903. }
  1904.  
  1905. void DetermineStructure()
  1906. {
  1907.     register Chain __far *chain;
  1908.     register Group __far *group;
  1909.     char buffer[40];
  1910.  
  1911.     if( !Database )
  1912.         return;
  1913.  
  1914.     if( InfoHBondCount<0 )
  1915.         CalcHydrogenBonds();
  1916.  
  1917.     if( InfoHelixCount>=0 )
  1918.         for( chain=Database->clist; chain; chain=chain->cnext )
  1919.             for( group=chain->glist; group; group=group->gnext )
  1920.                 group->flag &= ~(HelixFlag|SheetFlag|TurnFlag);
  1921.  
  1922.     InfoStrucSource = SourceCalc;
  1923.     InfoLadderCount = 0;
  1924.     InfoHelixCount = 0;
  1925.     InfoTurnCount = 0;
  1926.  
  1927.     if( InfoHBondCount )
  1928.     {   FindAlphaHelix(4,Helix4Flag);
  1929.         FindBetaSheets();
  1930.         FindAlphaHelix(3,Helix3Flag);
  1931.         FindAlphaHelix(5,Helix5Flag);
  1932.         FindTurnStructure();
  1933.     }
  1934.  
  1935.     if( CommandActive )
  1936.         WriteChar('\n');
  1937.     CommandActive=False;
  1938.  
  1939.     sprintf(buffer,"Number of Helices ... %d\n",InfoHelixCount);
  1940.     WriteString(buffer);
  1941.     sprintf(buffer,"Number of Ladders ... %d\n",InfoLadderCount);
  1942.     WriteString(buffer);
  1943.     sprintf(buffer,"Number of Turns ..... %d\n",InfoTurnCount);
  1944.     WriteString(buffer);
  1945. }
  1946.  
  1947.  
  1948. void RenumberMolecule( start )
  1949.     int start;
  1950. {
  1951.     register Chain __far *chain;
  1952.     register Group __far *group;
  1953.     register int hinit, minit;
  1954.     register int resno;
  1955.  
  1956.     if( !Database )
  1957.         return;
  1958.  
  1959.     hinit = minit = False;
  1960.     for( chain=Database->clist; chain; chain=chain->cnext )
  1961.     {   resno = start;
  1962.         for( group=chain->glist; group; group=group->gnext )
  1963.         {   if( group->alist->flag & HeteroFlag )
  1964.             {   if( hinit )
  1965.                 {   if( resno > MaxHetaRes )
  1966.                     {   MaxHetaRes = resno;
  1967.                     } else if( resno < MinHetaRes )
  1968.                         MinHetaRes = resno;
  1969.                 } else MinHetaRes = MaxHetaRes = resno;
  1970.                 hinit = True;
  1971.             } else
  1972.             {   if( minit )
  1973.                 {   if( resno > MaxMainRes )
  1974.                     {   MaxMainRes = resno;
  1975.                     } else if( resno < MinMainRes )
  1976.                         MinMainRes = resno;
  1977.                 } else MinMainRes = MaxMainRes = resno;
  1978.                 minit = True;
  1979.             }
  1980.             group->serno = resno++;
  1981.         }
  1982.     }
  1983. }
  1984.  
  1985.  
  1986. static void ReclaimAtoms( ptr )
  1987.     Atom __far *ptr;
  1988. {
  1989.     register Atom __far *temp;
  1990.  
  1991.     if( (temp = ptr) )
  1992.     {   while( temp->anext )
  1993.             temp=temp->anext;
  1994.         temp->anext = FreeAtom;
  1995.         FreeAtom = ptr;
  1996.     }
  1997. }
  1998.  
  1999. static void ResetDatabase()
  2000. {
  2001.     Database = CurMolecule = (void __far*)0;
  2002.     MainGroupCount = HetaGroupCount = 0;
  2003.     InfoChainCount = HetaAtomCount = 0;
  2004.     MainAtomCount = InfoBondCount = 0;  
  2005.     SelectCount = 0;
  2006.  
  2007.     InfoStrucSource = SourceNone;
  2008.     InfoSSBondCount = InfoHBondCount = -1;
  2009.     InfoHelixCount = InfoLadderCount = -1;
  2010.     InfoTurnCount = -1;
  2011.  
  2012.     CurGroup = (void __far*)0;
  2013.     CurChain = (void __far*)0;
  2014.     CurAtom = (void __far*)0;
  2015.  
  2016.     MinX = MinY = MinZ = 0;
  2017.     MaxX = MaxY = MaxZ = 0;
  2018.  
  2019.     MinMainTemp = MaxMainTemp = 0;
  2020.     MinHetaTemp = MaxHetaTemp = 0;
  2021.     MinMainRes = MaxMainRes = 0;
  2022.     MinHetaRes = MaxHetaRes = 0;
  2023.  
  2024.     *InfoMoleculeName = 0;    AdviseUpdate(AdvName);
  2025.     *InfoClassification = 0;    AdviseUpdate(AdvClass);
  2026.     *InfoIdentCode = 0;        AdviseUpdate(AdvIdent);
  2027.     *InfoSpaceGroup = 0;
  2028.     *InfoFileName = 0;
  2029.  
  2030.     VoxelsClean = False;
  2031.     HMinMaxFlag = False;
  2032.     MMinMaxFlag = False;
  2033.     ElemNo = MINELEM;
  2034.     ResNo = MINRES;
  2035.     MaskCount = 0;
  2036. }
  2037.  
  2038.  
  2039. void DestroyDatabase()
  2040. {
  2041.     register void __far *temp;
  2042.     register Group __far *gptr;
  2043.  
  2044.     if( Database )
  2045.     {   ReclaimHBonds( Database->slist );
  2046.         ReclaimHBonds( Database->hlist );
  2047.         ReclaimBonds( Database->blist );
  2048.  
  2049.         while( Database->clist )
  2050.         {   ReclaimBonds(Database->clist->blist);
  2051.             if( (gptr = Database->clist->glist) )
  2052.             {   ReclaimAtoms(gptr->alist);
  2053.                 while( gptr->gnext )
  2054.                 {   gptr = gptr->gnext;
  2055.                     ReclaimAtoms(gptr->alist);
  2056.                 }
  2057.                 gptr->gnext = FreeGroup;
  2058.                 FreeGroup = Database->clist->glist;
  2059.             }
  2060.             temp = Database->clist->cnext;
  2061.             Database->clist->cnext = FreeChain;
  2062.             FreeChain = Database->clist;
  2063.             Database->clist = temp;
  2064.         }
  2065.  
  2066.         FreeMolecule = Database;
  2067.         Database = (void __far*)0;
  2068.     }
  2069.     ResetDatabase();
  2070. }
  2071.  
  2072.  
  2073. void InitialiseDatabase()
  2074. {
  2075.     MinHBondDist = (Long)175*175;  MaxHBondDist = (Long)300*300;
  2076.     MinBondDist  = (Long)250*250;  MaxBondDist = (Long)475*475;
  2077.     AbsMaxBondDist = 475;
  2078.  
  2079.     FreeMolecule = (void __far*)0;
  2080.     FreeHBond = (void __far*)0;
  2081.     FreeChain = (void __far*)0;
  2082.     FreeGroup = (void __far*)0;
  2083.     FreeAtom = (void __far*)0;
  2084.     FreeBond = (void __far*)0;
  2085.  
  2086.     ResetDatabase();
  2087. }
  2088.  
  2089.